###  This script will try and produce Figure 1
### Figure one is three panels, WT DEM with site locs. BC05 diff and BC06 diff.  Needs continental Australia as well.  Scale bar and compass rose of some description.
### May need to be done individually and/or outside of R

### Out directory

out.dir = '/home1/99/jc152199/ChapterOne/'

### Load library

library('SDMTools')

### First read in the DEM and the model data file to get site locs

dem = read.asc('/home1/99/jc152199/MicroclimateStatisticalDownscale/250mASCII/STATIC/dem_WTplusbuffer_LatLong_WGS1984_250mres.asc')

model.data = read.csv('/home1/99/jc152199/brt/FINALOPTIMALMODELS/Model_Data_Plus_Preds.csv', header=T)
model.data2 = read.csv('//home1/99/jc152199/MicroclimateStatisticalDownscale/ToAnalyse/Model_Data.csv',header=T)

### Get a unique list of sites with xy

sites = unique(cbind(as.character(model.data$site),model.data$lat,model.data$long))
sites = data.frame(site=sites[,1],lat=sites[,2],long=sites[,3])
sites[,1] = as.character(sites[,1])
sites[,2] = as.numeric(as.character(sites[,2]))
sites[,3] = as.numeric(as.character(sites[,3]))

### Assign a column with site codes 1=CTBCC permanent 2=CTBCC intermittent 3=BoM

sites$pch = c(1,1,3,1,1,1,3,1,1,1,3,3,1,1,1,1,1,2,1,1,3,1,1,2,1,3,3,1,3,1,2,1,3,3,2,1,2,3,3,2,1,2,2,2,2,2,2,1,2,1,2,2,1)+15
sites[30,4]=2+15

### Assign an elevation code to each site, begin by merging DEM from model.data with sites

temp = unique(cbind(as.character(model.data$site),model.data$dem))
temp = data.frame(site=temp[,1], dem=temp[,2])
temp[,1] = as.character(temp[,1])
temp[,2] = as.character(temp[,2])

sites = merge(sites,temp,by.x=c('site'))

sites$dem = as.numeric(sites$dem)

sites$elevcode= 6

sites$elevcode[which(sites$dem<1250)]=5
sites$elevcode[which(sites$dem<1000)]=4
sites$elevcode[which(sites$dem<750)]=3
sites$elevcode[which(sites$dem<500)]=2
sites$elevcode[which(sites$dem<250)]=1

### Change DEM for site 32141 from NA to zero

sites[12,5]=0
sites[12,6]=1

##############################################################
### Remove site 32141

sites[nrow(sites),ncol(sites)]='#FFFFFF00'

###############################################################

### Create some color palettes

bnw.cols = colorRampPalette(c('black','white'))
greens = colorRampPalette(c('greenyellow','green4'))
blues = colorRampPalette(c('skyblue','slateblue4'))
reds = colorRampPalette(c('lightcoral','indianred4'))

### Add colors to plot to the sites dataframe

ttt=NULL

for (i in 16:18)

	{
	
	if(i==16) {cols=reds}
	if(i==17) {cols=greens}
	if(i==18) {cols=blues}
	
	tt = sites[which(sites$pch==i),]
	
	tt$col=NA
	
	tt$col = cols(6)[tt$elevcode]
	
	ttt = rbind(ttt,tt)
	
	}
	
sites=ttt

#### What follows is the script to produce an image of the WT DEM
#### Points will be overlaid for each climate sampling site
#### Characters are based on site type
#### Colors are based on elevation
#### Axes have been modified
#### No scale bar
#### There is a vertical gradient legend as well
	
	### Create an image of the DEM and overlay points

	png(paste(out.dir,'Figures/DEMwithSites.png',sep=''),units='cm',height=18,width=18,res=600)
	
	# Set plotting parameters
	
	par(oma=c(0,.5,.5,0), bg='#FFFFFF00')
	
	# Plot DEM

	image(dem, col=bnw.cols(100), bty='n', axes=F, xlab='', ylab='')
	
	### Plot axes
	
	axis(3, at=c(144.5,145.5,146.5,147.5), line=1)
	
	axis(2, at=c(-15.5,-16.5,-17.5,-18.5,-19.5), line=-5)
	
	### Use a for loop to plot a horizontal legend gradient

	### Establish x-y's for the corners

	pnts = data.frame(x=c(144.7,144.7,144.9,144.9),y=c(-17,-19,-19,-17))

	### Pick a number of colors from the bnw.cols palette

	dem.cols = c(bnw.cols(100))

	### Establish break points for gradient plotting

	yvals = seq(min(pnts[, 2]), max(pnts[, 2]), (max(pnts[, 2]) - min(pnts[, 2]))/(length(dem.cols) + 1))
		 
	for (inc in 1:length(dem.cols))

		{
		
		 polygon(x=pnts[,1],y=c(yvals[inc],yvals[inc+1],yvals[inc+1],yvals[inc]),col=dem.cols[inc],border=NA)
		
		}

	### Plot a hollow polygon with a thin black border around the DEM scale

	polygon(pnts,col='#FFFFFF00',lwd=.5)
	
	### Plot each site code with a different pch and color
		
	points(cbind(sites$long,sites$lat),col=sites$col,pch=sites$pch, cex=1)
	
	### Add some text labels
	
	text(144.8,-16.9,'1600',cex=.8)
	text(144.8,-19.1,'0',cex=.8)
	text(144.6,-18,'Elevation (m)',cex=.8,srt=90)
	
	# Flush output
		
	dev.off()
	
#### Now produce images of the difference between microCLIM and AWAPCLIM (only surfaces 5 and 6)

### List directories where surfaces are

mc.dir = '/home1/99/jc152199/MicroclimateStatisticalDownscale/250mASCII/microCLIM/'
awap.dir = '/home1/99/jc152199/MicroclimateStatisticalDownscale/250mASCII/AWAPCLIM/'

### List appropriate files in those directories

mc.files = list.files(mc.dir, full.names=T)[3:4]
awap.files = list.files(awap.dir, full.names=T)[3:4]

### Blank object 

tt = NULL

### Begin looping to create custom color palette

for (s in c('05','06'))

	{
	
	#### Select the surfaces to work with
	
	mc.file = mc.files[grep(s,mc.files)]
	awap.file = awap.files[grep(s,awap.files)]
	
	### Read in the surfaces
	
	mc = read.asc.gz(mc.file)
	awap = read.asc.gz(awap.file)
	
	### Calculate the difference
	
	dif = mc - awap
	
	#### Need to create a custom color palette that is grades from blue to white to red, white needs to correspond to a difference value of zero
	### When this loop is done, tt will be a vector of all possible difference values
	
	tt = c(tt,dif)
	
	}
	
### Find the range of tt

lims = range(tt,na.rm=T)
lims = c(lims[1]-.1,lims[2]+.1)
lims = round(lims,1)

### Create a custom color palette for plotting

### Start by establishing break points for colors

infrabreaks = seq(lims[1]-.1,-.1,.2) # Breaks below zero
ultrabreaks = seq(.1,lims[2]+.1,.2)	# Breaks above zero
allbreaks = c(infrabreaks,ultrabreaks)

#### Establish color palettes
#### NB Change white to a whitish shade
	
infracols = colorRampPalette(c('purple','blue','green','yellowgreen'))
ultracols = colorRampPalette(c('lemonchiffon','yellow','orange','red'))

### Get a list of colors based on the length of breaks from the color palettes

infracodes = infracols(length(infrabreaks)-1)
ultracodes = ultracols(length(ultrabreaks)-1)
allcodes = c(infracodes,'#FFFFFF',ultracodes) # This a list of colors one element longer than all breaks, with white lining up with zero

### Now loop through plotting

for (s in c('05','06'))

	{
	
	#### Select the surfaces to work with
	
	mc.file = mc.files[grep(s,mc.files)]
	awap.file = awap.files[grep(s,awap.files)]
	
	### Read in the surfaces
	
	mc = read.asc.gz(mc.file)
	awap = read.asc.gz(awap.file)
	
	### Calculate the difference
	
	dif = mc - awap
	
	### Open the .png device driver
	
	png(paste(out.dir,'Figures/FigureOne',s,'.png',sep=''),units='cm',height=18,width=18,res=600)
	
	# Set plotting parameters
	
	par(oma=c(0,.5,.5,0), bg='#FFFFFF00')
	
	if(s=='06')
	
	{
	
	# Plot difference

	image(dif, col=allcodes, breaks=allbreaks, bty='n', axes=F,xlab='', ylab='')
	
	}
	
	if(s=='05')
	
	{
	
	# Plot difference

	image(dif, col=allcodes, breaks=allbreaks, bty='n', axes=F,xlab='', ylab='')
	
	### Plot axes
	
	axis(3, at=c(144.5,145.5,146.5,147.5), line=1)
	
	axis(2, at=c(-15.5,-16.5,-17.5,-18.5,-19.5), line=-5)
	
	### Use a for loop to plot a horizontal legend gradient

	### Establish x-y's for the corners

	pnts = data.frame(x=c(144.7,144.7,144.9,144.9),y=c(-17,-19,-19,-17))

	### Establish break points for gradient plotting

	yvals = seq(min(pnts[, 2]), max(pnts[, 2]), (max(pnts[, 2]) - min(pnts[, 2]))/(length(allcodes)))
		 
	for (inc in 1:length(allcodes))

		{
		
		polygon(x=pnts[,1],y=c(yvals[inc],yvals[inc+1],yvals[inc+1],yvals[inc]),col=allcodes[inc],border=NA)
		
		}

	### Plot a hollow polygon with a thin black border around the DEM scale

	polygon(pnts,col='#FFFFFF00',lwd=.5)
	
	### Add some text labels
	
	text(144.8,-16.9,paste(round(max(lims,na.rm=T),1)),cex=.8)
	text(144.8,-19.1,paste(round(min(lims,na.rm=T),1)),cex=.8)
	text(144.6,-18,'Degrees C',cex=.8,srt=90)
	
	}
	
	### Flush output
	
	dev.off()
	
	}
	
# Done


	




	


